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ABSTRACT 

We analyze the stability of a relativistic double (forward/reverse) shock system 
j> ■ which forms when the fireball of a Gamma-Ray Burst (GRB) impacts on the 

surrounding medium. We find this shock system to be stable to linear global 
^vq , perturbations for either a uniform or a wind (r ) density profile of the ambient 

ON I medium. For the wind case, we calculate analytically the frequencies of the normal 

modes which could modulate the early short-term variability of GRB afterglows. We 
find that perturbations in the double shock system could induce oscillatory fluctuations 
^ ■ in the observed flux on short (down to seconds) time scales during the early phase of 

an afterglow. 

o 

CO 

Subject headings: gamma rays: bursts — ISM 

■ 1. Introduction 

Gamma-ray bursts (GRBs) and their afterglows are most naturally described by the 
relativistic "fireball" model (see e.g., Paczyfiski & Rhoads 1993; Katz 1994; Meszaros & Rees 1993, 
1997; Waxman 1997a, b; Sari, Piran, & Narayan 1998). In this model, a compact source releases 
a large amount of energy over a short time and produces a fireball which expands relativistically 
as a thin shell. When the shell encounters the circumburst medium, two shocks are formed - a 
forward shock which propagates into the circumburst medium and accelerates it, and a reverse 
shock which propagates into the relativistic shell and decelerates it (Meszaros & Rees 1992; Katz 
1994; Meszaros, Rees, &; Papathanassiou 1994; Sari & Piran 1995; Sari, Narayan, & Piran 1996). 
Later on, after a significant mass of circumburst medium is accumulated, the shell approaches a 
self-similar behaviour, as originally described by Blandford & McKee (1976), where there is only 
one forward shock propagating into the circumburst medium. The circumburst medium could be 
either the interstellar medium (ISM) or a progenitor wind. 
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The stability of the Blandford-McKee (1976) solution has been demonstrated recently by 
Gruzinov (2000). Here, we analyze the stability of a forward/reverse relativistic shock system. 
This double shock system exists during an important phase in the evolution of GRBs and its 
stability has observational consequences. In particular, oscillations or instabilities could translate 
to specific patterns of temporal variability in the lightcurves of GRB afterglows. 

In our linear perturbation analysis we generalize the "thin shell" method first introduced by 
Vishniac (1983) in the non-relativistic regime. This method simplifies the equations describing 
the stability of a spherical shock when the wavelength of the perturbation is much larger than the 
thickness of the shocked shell. In our relativistic treatment we focus on global perturbations for 
which the wavelength is much larger than the thickness of the forward/reverse shock system. We 
consider the regime of GRB parameters where the reverse shock is relativistic (although in reality 
it may also be non-relativistic). In §2 we derive the perturbation equations for the forward/reverse 
shock system. In §3 we show the analytical results for the wind case and the numerical results for 
both the wind and ISM cases. Finally, in §4 we summarize our main conclusions. 



2. Linear Perturbation Equations 

The interaction between a relativistically expanding shell and the circumburst medium results 
in a double shock system, as shown in Figure 1. The system includes four distinct regions: the 
circumburst medium (region 1) and the shocked circumburst medium (region 2) are separated by 
the forward shock (shock 1), while the shocked material in the shell (region 3) and the unshocked 
material in the shell (region 4) are separated by the reverse shock (shock 2). Region 2 and 3 
are separated by a contact discontinuity. Our analysis is done in the observer frame where the 
circumburst medium is at rest. We use a spherical coordinate system whose origin is located at 
the center of the explosion. The radii of the contact discontinuity, shock 1 and shock 2 are denoted 
by Ro, R\ and i?2 respectively. We refer to the combination of region 2 and 3 as the layer whose 
stability we consider. Similarly to Vishniac (1983) we make the thin shell approximation, i.e. 
assume that 

p < HRl ~ R2) « 1, (1) 

where k is the wavenumber of the perturbations. Note that although shock 1 is relativistic, shock 
2 could be either relativistic or nonrelativistic. In this paper, we only consider the situation where 
shock 2 is relativistic. 

The equations of motion for an ideal relativistic fluid are 

_( 7/0 ) +V .( 7/ 9u)=0, (2) 



7 2 



<9 U / ' 
^ + (u.V)u 
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where p, u, e, p and 7 are the fluid density, velocity, energy density, pressure and Lorentz factor, 
respectively. We define the surface density a, bulk radial velocity V r and average tangential 
velocity Vy of region 2 as follows: 

a(9, ( f>)=R 2 [ R \pr 2 dr, (4) 

V r (0, <f>) = {<tRIY 1 1P u r r 2 dr, (5) 
JRo 

'■Hi 
IRo 

The time evolution of these variables can be obtained by integrating equations @ and (||) across 
region 2, using the boundary conditions at shock 1 and at the contact discontinuity and neglecting 
terms of higher order in (R\ — R^/Rq an d k(R\ — i^)- We get 

d t a = -2-j^cr + pic - cjV t • V r , (7) 
8a(Ro) = -2a-HR»)P,c + + ^VrV T • V Tl (8) 



V T (0, 0) = (cT^)" 1 7 pu T r 2 dr. (6) 



„ _o3/4 o l/2 P 3/4 (^0)Pl /4 
Pl 7 3/ 2 (i?0)c 3 / 2 



VtRq — c 1 / icVy 



-^V T - a" 1 — 1_ V T <r + -^—V^o), (9) 
#0 37 2 (i?o) 37^(^0) 

where pi is the density of the unshocked circumburst medium just in front of shock 1, 7 (i?o) 
and p(Ro) are the Lorentz factor and pressure at the contact discontinuity, and Vy denotes the 
tangential derivatives. In deriving the above equations we also assumed that the radial velocities 
are dominated by the bulk motion of regions 2 and 3 (denoted thereafter as the "shock layer" ) , so 
that Rq ~ V r . The full derivation of the above equations is given in the Appendix. 

From equations (0) and (|8|) we obtain the unperturbed equations: 

fta<°> = -2^°) + pic, (10) 
R 

q3/4 3/4 1/4 

a 7c = - 2 (a(«»)-w+ w M 0) r lP ^r^ (id 

where u c = i?g°^ = V^ ^ and 7c = 7^(-Ro) = — v'^/c 2 are the velocity and Lorentz factor of 

the unperturbed contact discontinuity, p c = p'°'(i2o) is the unperturbed pressure at the contact 
discontinuity, and we use a superscript (0) to denote unperturbed values. 
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For the shocked shell in region 3, we define the surface density, bulk radial velocity and 
average tangential velocity to be: 



JRo. 



V r3 (e,0) = (a 3 R 2 o 



2\-l 



Ro 

R-2 

Ro 



^pu r r dr, 



\ T3 (0,<t>) = (<T 3 R 2 )- 1 7PUT 

Jr 2 



r 2 dr, 



(12) 

(13) 
(14) 



where a subscript 3 denotes quantities in region 3. The time derivatives of the above variables can 
be derived similarly to those in region 2, 



d t a 3 = -2^-o 3 + 2 A P4C - a 3 V T • V T3 , 

iP 3/4 («o)P4 /4 7 1/2 («o) , i(Ro) 



d a (R )=a^^- Pi c-3 3 / 4 a- 3 
l(Ro) 



(15) 



+ ^_^y r V r -V T3 , (16) 
2c 



d t V T3 



2.38/4 P 3/ HRo)pT 



74 



7 4 1/ V/2( j R ) c 3/2 272(^0) 



Pi 



VtRo 



-(To 



-1 74 



l 2 (Ro 

o2 



v r 

-p 4 cV T3 - — V T3 



4 7 2 (^o) V 



V T3 



rr^ 1 - o e „ , Vro- 3 + - q/n Vr7(-R ), 



3 7 2 (#o) 



37 3 (^o 



(17) 



where p^ and 74 are the density and Lorentz factor of the unshocked shell (region 4) just in front of 
shock 2. We have also used the relation V r3 ~ V r , as appropriate in the thin shell approximation. 



The unperturbed equations for region 3 are 



a (°) 



o Vc (°) , 74 



^7c = (4 0) )- 1 ^P4C-3 3 / 4 ( 

7c 



(7 



3/4 1/4 1/2 

(0)^-1 Pc_ Pa lc 

1/2 1/2 

7 4 c-V<* 



(18) 



(19) 



Equations (Q), @, (|9|), (|l5|), (16) and (17) make a complete set of equations for the evolution 
of the forward/reverse shock system. The perturbation equations depend on the density profile of 
the circumburst medium. In the following subsections, we shall derive the perturbation equations 
for a uniform medium (such as the ISM) and for a progenitor wind. 
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2.1. Uniform Medium (ISM) 



For a uniform circumburst medium, p\ = const. We define the perturbation variables 
5 = a/aM -1,S 3 = o- 3 /4 0) - 1, Aft = Rq — R^ and A p = p(R ) / 'p (0) (ft ) - 1. From equations 
(|7|), (||) and (^), we obtain the following perturbation equat 



tions: 



8,8 



R 



(0) 



d t AR + 2 V " AR - (aW)- 1 Pl c5 - V r • V r , 
(4 0) ) 2 



(20) 



dfAR 



2( ff (°))- 



+ 



4(^°)rVic 



Q3/4 3/4 1/4 1/2 

|_ (t7 (0))-i ^ Pl c 
7 ■ 3 3 / 4 , 



+ 



3 7 / 4 
2 5 / 2 



23/2 
3/4 1/4 1/2 



a 



7/2 

7c 

3/4 1/4 
(0)j-lPc Pl 



3/2.1/2 



fr (0) )- l Pc V A P + ^V T -V T , 

7c 



7c c 
27 



(21) 



9 t V T = \{a^)- l c 2 



Pi 



3/4 1/4 
g3/4 q1/2 ^ c Pl 

7c 3/2 C 3 / 2 



V T AR- (a (0) )"VicV T 



ft 



^V T -^V T * + fv T (ftAl2). 



(22) 



Assuming that 74 is a constant and there is no shell spreading, we get P4 oc R 2 . Using this 
scaling we derive the following perturbation equations from equations (p!q), (|l~6|) and (|i~7|): 



■t0 3 



-2(4° ) )~ 1 74P4^Ai? + 



(0)\2 



, (0)s-i74 x 
- (°3 ) -5P4C63 

7c 



L (ftT) 



2(4 u; ) 



(0)^-1 74 P4C 



2 p(°) 



7£ R 



J 



Ai? 



(23) 



9 2 Aft 



-(4° ) )- 1 ^P4C- + ^^CT3 

7c 



.2 - o3/4/ (OK- W Pi C l/ 
+ 6 V°3 J 1/2 5/2 

7 4 7c 



+ 



+ 



i^ (0) )- l74 ^c 1 5 ' 33/4 ( J (0) r 1 p ' /4 ^ /4 

4 ^3 J ^2^ + 9 13 ' 1/2 1/2 I/a 
7c z 7 4 ' 7c' C 1 /^ 

, J0K-1 74 P4C 2 3 3 / 4 (Q) ! Pc /4 ^ /4 c 1/2 1 

- !rT 3 1 riTTir • - ^ _ifJ :; i 



72 R 







1/2 5/2 
7 4 7c' 



R 



d f AR 



AR 



(0) 



3 7 /4 ... 3/4 1/4 1/2 



1/2 5/2 

7 4 7c' 



27 2 



V 



T3, 



(24) 
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dt v T3 = l(4 0) )-v 



3/4 1/4 

^V'W o.:VVi_L y " 

L / J /'• 

i /<W 



9 q3/4 p 4 ' _ 74 

1/2 1/2 3/2 2>y 2p4 
74 7c c J /^ Z 7 C _ 



7c 



V T3 - -r-5 ( — J V ^3 - ^V T <5 3 + ^V T (3 4 Ai?). (25) 



^(o) 4 7c 2 v P4 ; 3 7c 2 



In total, we have six perturbation equations p0j)-(p5[) in six variables: 5, 8s, AR, Vy, V^3 
and A p . In order to solve these equations, we first need to find the unperturbed values 
7 C and p c from equations (JlO|) , (|TT|) , ( |18D and (|1S|). The time dependence of 7 C has been derived by 
Sari & Piran (1995) and Sari et al. (1996). When both the forward shock and the reverse shock 
are ultrarelativistic and strong, 

7c oc 7 4 1/2 / 1/4 , (26) 
where / = Pi/ pi- For p\ = const, 74 = const and p± oc R~ 2 , we get 



7c oc 



R- 1 ' 2 oc t- 1 ' 2 . (27) 



With v c m c and R^ « ct, equation (|10|) yields 



^ (0) « -pict. (28) 



For P4 oc t 2 and j c oc t 1 / 2 , equation (18) gives 



o"^ « ~--|y04ct = const. (29) 
^ 7c 



By substituting equation fl28j) into equation (|ll|) we find 



1 1 4/3 4 

2 2 1 1 07 2 2^*2 2 /nn-, 

7cPic z = 1.1877cPic z < -%pic', (30) 



^ C 3 7 /3. 2 2/3 ' c ^~ ' ^ 3 

where (4/3)7 2 /9ic 2 is the pressure just behind the forward shock. Similarly, by substituting 



equation ( |29| ) into equation (|1^) we get 

5 4 / 3 74 2 n ^nT4 2^4-22 . , 

where 73 ~ 74/(27,;) is the Lorentz factor of the shocked shell (region 3) with respect to the 
unshocked shell (region 4) and (4/3)7 3 2 p4C 2 is the pressure just behind the reverse shock. The 
pressure difference between the two sides of the layer causes it to decelerate. By combining 
equations ( |30[ ) and ( |3l| ) we find 

7c« 0.7847 4 1/2 ( P4 /pi) 1/4 . (32) 
Substitution of the values of and p c into the perturbation equations (pc|)-(|25|) yields, 

dt5 = ~atAR+-^AR-*8-VT-VT, (33) 
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a ' AR = lit 6 - fl 9 ' AB + f + if T ■ Vt - (34) 

8,Vt = -5£v T Aii - jV T - f^V T S + §V T (ftA,R), (35) 
3 t t 07^ 3 

d t 5 3 = ~^d t AR - ^AR - ~£ 3 - V T ■ V T3 , (36) 

= ~h - l~ t O t AR ~ H-^Afl - H -E.A, + ^V T • V T3 , (37) 

QtVra = —VtAR - ^V T3 - V T 5 3 + ^V T {d t AR). (38) 
3 1 t Sjz 3 

Combining equations ( |3^ ) and ( |37| ) and eliminating A p , we get 

1 32 7 2t 32 Ht 6 32 r 64 7^2 



+ |^V T .V T + | : |v T .V T3 . (39) 



Next we expand the spatial dependence of the perturbation variables in spherical harmonics. 
We choose to normalize these variables so as to make them dimensionless and of a similar 
magnitude through the definitions 

AR = J2 Ai*(Z, m, t)[R^ hl\Y lm {9, </>), (40) 
6 = ^28(l,m,t)Y lm (0,<i>), (41) 

l,m 

= m,t)lW(M), (42) 

V T =J2V T (l,rnMcli 0) hc] VTYl f^ , (43) 



V T3 = ^V T ^ m Mc^Vlc] VTYl f^ . (44) 



Equations (|33|), (pq), (pq), (pq) and ( p9| ) can be rewritten as 

d t AR = F, (45) 



= -- 2 F- \AR - h + ^-V T , (46) 
7c 7c* * Id 



1 11 1 Q 

W = - — Ai? - - -F T , (47) 

37 c j c t 3j c t 2t 
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d t 5 3 = -AF - -AR - ~<5 3 + ^P-Vts, (48) 

1 21 17 

d t V T3 = - — F + — AR - —5 3 - -Vt3, (49) 
37 c 7 c t 3j c t 2t 

239 ^ 143 5 . 11 . 5(Z + l) Tr 11 (Z + 11 

d t F = F kAR-\ ~5-\ k5 3 t^Vt ^VW (50) 

32t 16t 2 32t 2 32t 2 A 32 7c t 2 32 7c t 2 6 y 1 

The perturbation variables in the above six equations are dimensionless and only functions of 
time. We have added the variable F so that all the equations will have the form of first order 
differential equations. 



2.2. Wind Medium 



If the circumburst medium is a progenitor wind, p\ oc R 2 . Accordingly, the perturbation 
equations (|20|) and (|2l|) need to be changed to 



d,S 



R 



(0) 



d t AR + 



'(4 0) ) 



2(a 



(Q)x-l P\c 



R 



(o) 



AR - (ct (0) )-Vic5 - V T • V T , 



(51) 



2(a^) 



(0)\-lPl c 



3 3/4 



7, 



2 2 V2 



cr 



7/2 
7c 



+ 



+ 



7 o3/4 r, 3 / 4 ^ 1 ^ 



4(fj 



(O)x-l Pl c 



3 3/4 



7 2d(0) 2 3 / 2 

7c -"-o 



, ^ D 3/ V /4 r-V 2 
(0)l-l Pc Pi c 

7/2 (0) 

7c -ftQ 



AR 



Q7/4 3/4 1/4 1/2 

+ 2 5/2^ J 7 J/2 27 2 

Equations (p2^)-(|25|) remain the same as in the uniform medium case. 



(52) 



Since p\ oc R 2 and p^oc R 2 in the wind case, equation (^) implies that 7 C is constant over 
time. Equations (10), (11), (18) and (19) then yield the unperturbed parameters: 



a® « pxct, 
7c 

4 2 2 1-74 2 

« ,7cPl c ~ o— P4C , 

7c 

1 1/2 fp^ 1/4 



(53) 
(54) 

(55) 

(56) 
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Substitution of equations (p3|)-(pq) into the perturbation equations (51), (52) and (p2|)-(25), 
yields 

(57) 
(58) 

(59) 

(60) 
(61) 

(62) 



dfAR 



d t 5 = -—d t AR -S-V T - V t , 
ct t 

3 , . . 3 . 3c. c , _ ,. 

iatAR+ _ AR+ ._. Ar+ _ VT . VT , 

-ix7 T AR - ly T - f~V T S + ^Vr(9,Afl), 
t t djg 6 



d 2 t AR 



31 



-2^d t AR - 1*3 - V T • V T3 , 
ct t 



2i 9 ' AR -lw AR -Wt A ' + Wf T ' VT3 - 
1 c. 



V T Ai? - -V T3 - — j V r <$3 + -V T {d t AR). 
2 1 t S"j£ 6 



9fV T3 

Combining equations ([58|) and fl6l|) , and eliminating A p , we get 

d 2 Ai? = -H^AB + V r • V T + -^V T • V r3 . 

t 67^ 37^ 



(63) 



Using the same normalized perturbation variables as defined in equations (|Ic|)-(|44]), we can 
rewrite equations (|57|) , (|59|), (|60|) , (|62]) and (p3) as 







d t Ai? = F, 


(64) 




-V- 

7c 


7c* t let 


(65) 




= -Lf 

37 c 


21 AR 1 5 2 V T , 
37ct 37 c t t 


(66) 




= -2F- 


l AR .l S3 + C^v T3 , 

t t 7ci 


(67) 


d t V T3 


Z 

= F- 

37 c 


5/ I 2 

^ Ai? <5 3 V T 3, 

67 c t 37ct t 


(68) 


d t F = - 


If- ' 
t t 2 


a*- ( ' + 1 2 V t - (/ + 1 2 V t , 

6 7c t 2 3 7c t 2 


(69) 



These six first order differential equations are the final perturbation equations for the wind case. 
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3. Solutions of the Perturbations Equations 



For the wind case, j c = const, and we can solve the perturbation equations analytically. If we 
define F = F' /t, then equations (|64|)-(|69|) can be rewritten in a matrix form: 



( F ' ) 




( 


-3 


-2 





67c 





(1+1) 

37c 


AR 






1 

















5 


1 




2 


2 


-1 


(i+i) 








V T 


~~ t 




I 

37c 


21 

37c 


l_ 

37c 


-2 














-2 


-2 








-1 


(1+1) 

7c 


\ V T 3 J 




\ 


i 

37c 


_5l_ 

67c 








i_ 

37c 


-2 









AR 




5 




V T 








\ VT3 J 



(70) 



In matrix notation, the above equation is 

d t y = • y, 

where y is a vector, and A is a 6 x 6 time-independent matrix. 

The matrix A can be diagonalized through the transformation 

diag(Ai • • • A 6 ) = D, 



(71) 



X" 1 • A • X 



(72) 



where Ai, • • •, A6 are the eigenvalues of the matrix A, and X is the matrix formed by columns 
from the eigenvectors (i.e. the £;-th column of X is the eigenvector corresponding to the eigenvalue 
Afc). Equation ( f7l"l) can then be transformed to 



d t y = -(X • D • X- 1 ) • y =► ^(X" 1 • y) 

By defining a new vector y' = X -1 • y, we get the equation 

I. 



-D (X 



•yj 



d t y' = -D y' 



t 



(73) 



(74) 



which has six components 



1 



d tVk = jhVk, 



fc = l,---,6. (75) 
Since A^ can be a complex number, we write A^ = a& + ibk- The solution to equation ( |75| ) is then 

c k t ak+ibk = c k t ak e ibklnt , (76) 



y'k 



where c k is a constant dictated by the initial conditions. Each y' k defines a mode of the shock 
system. There are six modes in total corresponding to six eigenvalues of the matrix A. The real 
part of each eigenvalue dictates the overall temporal behavior of the mode, while the imaginary 
part determines its oscillation frequency. The vector y can be derived from the relation 



y = X • y'. 



(77) 
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Hence, each component of the vector y is a linear combination of the six different modes. The 
vector c whose components are cjt, can be obtained from the initial conditions 

c = y / (t = l)=X- 1 -y(t = l) 

namely 



(78) 



( 01 ) 




/ F'(t = l) \ 


C2 




AR(t = 1) 


C3 


= x- x - 


S(t = 1) 


C 4 


V T (t = 1) 


C5 




= 1) 


V C 6 ) 




V V T 3(t = 1)) 



(79) 



The eigenvalues of the matrix A can be calculated for different values of "y c and I. Figure 2 
shows the real and imaginary parts of the six eigenvalues as functions of I for 7 C = 500. Each row 
in the figure contains two panels, corresponding to the real and imaginary parts of a particular 
eigenvalue. The results show that for I < 320, all the six eigenvalues are real numbers, and so 
there are no oscillations. For 320 < / < 432, two eigenvalues are complex numbers and they are 
a pair of complex conjugates, implying that two modes are oscillating with the same frequency. 
For / > 432, there are two pairs of complex conjugates. The transition from real eigenvalues to 
complex eigenvalues occurs when I ~ j c , as expected from the fact that oscillations are possible 
only when causality allows communication across the scale of a wavelength for modes with / ^> 7 C . 

For the thin shell approximation to be valid, we require that the wavelength of the 
perturbation be much larger than the thickness of the forward/reverse shock system in the shock 
frame. The thickness of the shock system is ^ 2i?o/C7c m the observer frame, and thus <^ 2i?o/C7c 
in the shock frame. Here C is a constant that ranges between ~ 4 and ~ 12 for the wind and 
ISM profiles, respectively. The wavelength of the perturbation is ~ 2itRq/1 in the shock frame. 
Therefore, we enforce an upper limit on / of ~ 107 c . Figure 2 shows that for all values of I, the 
real parts of the six eigenvalues are < —1. This implies that all modes are decaying faster or 
proportional to t . Since each perturbation variable is a linear combination of the six modes, we 
conclude that all perturbation variables should also decay faster than or equal to t . Thus the 
system is stable. Note that for large 7 C we expect the results to depend only on l/j c [see equation 
(|70|)1, and so our particular choice of 7 C = 500 can be scaled appropriately to other values of j c . 



For / 3> 7 C the eigenvalues admit the following analytical solutions, 
Ai = -1, A 2 



3 11 

-2-vSTc'- 



3 1 I 
— I — i, 

2 ^/3 7c ' 



. 19 13 1 I . 13 1 I . 
Aa — , A5 — 1 — — i, Afi — — — i, 

9 9 V27c 9 V27c 



while in the limit of I <C 7 C , 

Ai 



1, A 2 



-2 + 



IP 

3 7c 2 



1P_ 

37 C 2 ' 



(80) 
(81) 

(82) 
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1 I 2 11, 11 

A 4 = -l A 5 = -2 =— , A 6 = -2 + — =— . 83) 

2 7 C 2 3V2 7c 3V2 7c 

We have calculated the corresponding eigenvectors numerically as shown in Figure 3 (for 
7 C = 500). Each row in the figure contains two panels which show the real and imaginary parts of 
one of the six components of the eigenvectors as functions of I. Different line types correspond to 
the six different eigenvectors. The complex eigenvectors are all scaled to have a unit magnitude. 
Since each eigenvector corresponds to a mode, the relative values of the six components of the 
eigenvector measure the physical significance of perturbations in different physical parameters for 
that mode. For example, the mode corresponding to the eigenvalue Ai = — 1 has the temporal 
behavior of t ; the eigenvector for this mode is (-0.275, 0.275, -0.824, 0, 0.412, 0), implying that 
this mode does not involve Vr and Vt3 perturbations. Also note that this mode does not depend 
on / while all other modes change with I. 



Equations (p4|)-(|69|) can also be solved numerically. We normalize all initial values of the 
perturbation variables to unity. The temporal interval of the calculation is from t = 1 to 30, 
and 7 C is chosen to be 500. The time t = 1 marks the beginning of the double shock system. 
After the initial explosion, the fireball expands and accelerates to a relativistic speed. A cold 
shell is formed after the thermal energy of the fireball is converted to a bulk kinetic energy. The 
forward/reverse shock system is formed at radius Rt=i where the circumburst medium starts to 
decelerate significantly the relativistic shell. This radius can be approximated as the radius at 
which the mass of the circumburst medium swept by the shell is comparable to Mo/77, where Mq 
is the mass of the initial baryonic load of the fireball, and r/ is the initial thermal Lorentz factor of 
the fireball. If the fireball has a total equivalent isotropic energy E, we have Mo = E/r]c 2 . Thus 
t = 1 corresponds to a time Tt=\ ~ (1 + z)Rt=\j1^ 2 c after the GRB trigger in the observer frame, 
where z is the cosmological redshift of the source. Figure 4 shows our results. The six panels show 
the time evolution of A.R, d t AR, 5, Vr, 63 and Vt3- We show results for four different I values, 
namely I = 5, I = 50, I = 500 and I = 5 X 10 3 . These plots indicate that all perturbation variables 
decay quickly with time. For small values of I (e.g., I = 5 and I = 50) there are no oscillations. For 
large values of I (e.g., I = 5 X 10 3 ) the oscillations exist but damp away quickly. These results are 
consistent with our analytical derivations. For I = 500, the oscillations start to appear although 
they are not apparent in the plot because of their low frequency. For I 3> j c , we can calculate 
the frequencies of the oscillations using the eigenvalues listed in equations ( |80| ) and (^). The 
oscillations have the form exp(i^^- hit) or exp(i-^^ hit), so that the oscillation period increases 
with time. For I = IO70 the shortest period is ~ 2.4, corresponding to 2.4Tj = i in the observer 
frame. Hence, the oscillations could produce fluctuations in the observed flux on time scales as 
short as a few times the starting time of the double shock system Tf = \. 



In the ISM case, 7 C is not constant and so we can not write equations (f45|)-(|50|) in a matrix 
form that admits an analytic solution. Instead, we had to solve these equations numerically. In 
order to test the validity of the perturbation equations and the numerical code, we compared 
the numerical results for a spherical perturbation with I = to the analytic solution derived 
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by perturbing directly the radial equations of motion, and found an excellent agreement 
between the two calculations. Similarly to the wind case, we normalized all initial values of 
the perturbation variables to unity, and chose an initial 7 C = 500. Our numerical results are 
shown in Figure 5 and resemble qualitatively the wind case. Overall, the perturbations decay 
rapidly with time and oscillations appear only for large values of I. Similarly to the wind 
case, for I = Wj c the shortest period of the oscillations is ~ 2, corresponding to ~ 2Tt=\ in 
the observer frame. If the equivalent isotropic energy of the fireball is E = 10 52 ergs, the 
initial thermal Lorentz factor of the fireball is r] = 10 3 , the number density of the ISM is 
n = 1 cm -3 , then we have Rt=i = (3£ , /47rnm p ry 2 c 2 ) 1//3 ~ 1.2 x 10 16 cm. This corresponds to 
T t= \ = (1 + z)Rt=i/2^c ~ 0.8(1 + z) sec. Thus, the time scales of the fluctuations in the observed 
flux could be as short as ~ 2(1 + z) sec. 



4. Discussion 

We have solved the perturbation equations describing the double (forward/reverse) shock 
system which forms during the impact of a highly relativistic fireball on a surrounding medium. 
For both a uniform and a wind (1/r 2 ) density profile of the ambient medium, we have found 
the shock system to be stable to global perturbations. We therefore do not expect the shock to 
fragment. Our results are limited to relativistic reverse shocks, and appear to differ qualitatively 
from previous results in the non-relativistic regime (Vishniac 1983). 

Our results apply also to collimated outflows as long as the double shock system is formed at 
a time when the Lorentz factor of the outflow is larger than the collimation angle. 

We derived the frequencies of the normal modes which could modulate the short-term 
variability at the early phase of GRB afterglows. The results imply that perturbations in the 
double shock system could produce fluctuations in the observed flux on time scales as short as a 
few seconds for j c ~ 500 in the ISM case. These short term fluctuations could be supplemented 
by variability on much longer time scales due to density inhomogeneities in the ISM; such 
inhomogeneities can lead to variability on time scales of tens of minutes in the optical band and 
days in the radio (Wang & Loeb 2000). 

This work was supported in part by grants from the Israel-US BSF (BSF-9800343), NSF 
(AST-9900877; AST-0071019), and NASA (NAG5-7039; NAG5-7768). 
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APPENDIX 

Here we provide full details for the derivation of equations in §2. We start by listing 

the equations of motion for a relativistic fluid in spherical coordinates. The continuity equation 
reads 

did Id Id 

oi(7P) + —Tr( r2 lpur) + ——^(smejpue) + ———Upu^) = 0, (1) 
dt r £ dr r sin da r sin o<p 

and the three components of the momentum equation are 

7 2 du r du r ldu r 1 du r 1 „ 2 dp u r dp 

ot or r 90 v rsm6 d(j) r v dr c 1 dt 

l 2 , . „du e du e . I du e _ldu e 1 2 1 dp 

c 2 ot <9r r w rsm# dq> r v r do c 2 dt 

T 2 / . . ^ , 1<9tu 1 dut I. 1 dp 114 dp . . 

— (e+p)[^+n r — +u e -^+u^^— ^+-(« r ^+cot6lu 9 u )J + — — — + — — = 0, (4) 
c z at ar r « r sin a0 r r sm a<p cr ar 

where /9, e, p and 7 are the fluid density, energy density, pressure and Lorentz factor respectively, 

u r , ue and u<f, are the three components of the fluid velocity. 

For the forward/reverse shock system under consideration (see Figure 1), we define the 
following shell-averaged variables for region 2: 

a(6,0)=R Q 2 f R \pr 2 dr, (5) 
JRo 



V r [6,<f>) = (oRlY 1 f 1 1P u r r 2 dr, (6) 
JRo 

V r (0, 0) = (aR 2 )" 1 f Rl 1P u T r 2 dr, (7) 
JRo 

where uy is the tangential velocity vector. 

Since the shocked material is relativistic, we adopt the relativistic equation of state, p = e/3, 
in region 2. Equation (|l]) yields the evolution of the surface density 

= -2^a + (j^j 1 {R 1 )p{R l )[R 1 - Ur(Rx)] + j(R )p(R )[u r (Ro) - Ro] 

-R 2 [ Rl [V-( 7 pu T )]r 2 dr. (8) 
JRo 

Since r = Rq defines the contact discontinuity between the shocked shell and the shocked 
circumburst medium and there is no mass flow across the contact discontinuity, we get 
u r (Ro) = Ro. Because the forward shock is a strong relativistic shock, we have the following shock 
jump conditions at shock 1: 

l 2 (Ri) = 7 2 i/2, (9) 
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p(Ri)/ Pl = 47(i?i), (10) 

where j(Ri) and p(Ri) are the Lorentz factor and density of the fluid just behind the shock front, 
7 S 1 is the Lorentz factor of the shock front and p\ is the density of the unshocked circumburst 
medium just in front of the shock front. In the highly relativistic regime, 



1 

Hi 



Ri~c(l-— , (11) 



0(1-5^). (12) 

From equations (p|)-(|l2"|), we obtain 

7 (i?i)p(i?i)[i?a - Ur(Ri)] « pic. (13) 
Thus, equation (|8|) can be rewritten as 



To linear order, the last integration term in the above equation can be approximated by 

Ro 2 [V • (7Pu T )]r 2 dr = -ctVt ■ V T + R 2 [** [V • (jpu T )) - l) r 2 dr, (15) 
Jro JRo \Ro J 



where the operator Vt = (9/Ro)(d/d6) + ((f) /Rq)(8/ d<j>) acts as follows on a scalar \£ and a vector 
f: 

1 - 1 - 
V T tf = — ^0 + ^U, (16) 

ito sm 6* aw i?o sin a<p 



Note that the second term on the right-hand side of equation (|lq) is of higher order in — Rq)/Rq 
than the preceding term, and hence can be ignored in the thin shell approximation. Thus equation 
(14) can be rewritten as 

d t a = -2^a+(^J pic-cjVt-Vt. (18) 
Similarly to the above derivation, we obtain for the bulk radial velocity 



d t v r (e,</>) = -v Fr_2 iV r + a (i^J ^( R i)p( R i) u r(Ri)[Ri-MRi)} 



+a- 1 j(Ro)p(Ro)u r (Ro)[u r (R ) - R ] - {aRlT 1 [ V • {7 pu T )]u r r 2 dr 

JRo 

-<"*>- 1 £ £ (£ + s t) r2dr - £ w(ut • v "* 2 * 

+ (a^)" 1 / jpu^rdr. (19) 
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The last two terms on the right-hand-side of the above equation are nonlinear. By substituting 
equations ( |l3| ) and (|14]) into the above equation, and keeping terms to the linear order we get 



d t V r (9, 



a 



+ 



-1 



R\ \ 2 ( , , 
—j p 1 c[u r {Rx J 



V, 



Ro 
«i pc 2 



[V • ( 7P u T )]r 2 dr - (aRlY 



Ri 



Ro 



u r[V • (7 pux)]r 2 dr 



dp u r dp 
Ro ^1P V^r c 2 dt 



r dr. 



(20) 



In order to evaluate the integral in the last term of the above equation, we need the relation 
between p and p inside region 2. Since entropy is conserved in this region, 



d_ 

dt 



P 

4/3 



0, 



(21) 



implying that p/p 4 / 3 remains constant for a given fluid particle. Hence, a fluid layer which is at a 
distance x(t) from the contact discontinuity (r = Rq) inside region 2, maintains a constant p/p^ 3 
over time and its value is decided by the Lorentz factor of shock 1 at the time when this layer 
first crossed shock 1. However, at a particular time, different layers across region 2 have different 
values of p//? 4//3 . Assuming that region 2 is decelerating with 7 oc r -1 / 2 (for the uniform ISM 
case), we get 



p(x) 



1/3 



E>V 2 
laRa 



-> 2/3 



(8riR a x + R 2 )^ 



(22) 



where j a and R a are the Lorentz factor and radius of region 2 at the initial time. Apparently, 
the dependence of p(x) / p^ 3 (x) on x is very weak, and so within the context of the thin shell 
approximation we simply assume that p/p 4 / 3 is constant across region 2 at any given time. This 
assumption is indeed satisfied in the numerical simulations performed by Kobayashi, Piran, & 
Sari (1999). In equation (|^), the term 8j 2 R a x can be at most comparable to R 2 (this happens 
in the very last stage of the evolution when the reverse shock crosses the shell), and so the error 
introduced by our approximation is small. For the wind case, 7 ~ const, and we can also treat 
pj p 4 / 3 as a constant across region 2. 



We may now calculate the integral 



Ri 



pc 2 dp 2 

r dr 



r Ajp dr 



Ri 



C d P 2 , 

— — r dr 

Ro 37 dr 



37(^o) 



Hi 

Ro dr 



d P 2 W 
— — r dr 



37(^0 



piR^R 2 - p(R )R 2 



Ri 



Ro 



2rpdr 



(23) 



In the thin shell approximation, all radial velocities are dominated by the overall radial motion of 
the shock layer. Hence, 7 was treated as a constant and was taken out of the integral. For the 
uniform ISM case, the density difference between the two edges of region 2 is not small; hence, 
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the last term inside the square brackets is of order (R\ — Rq)/Rq times the difference between the 
previous two terms and so it can be neglected. For the wind case, there is no pressure gradient 
across region 2, and so the above integral vanishes. However, even in this case we can ignore the 
last term and keep only the first two terms, because later on we will replace R± with Rq, and so 
the difference between the first two terms will vanish. 



Another relevant integral is 



Rl pc 2 u r dp o 

r dr 



r 47p c 2 dt 



fl l C 2 U r dP r 2 dr _ V r f Rl d P r 2 dr 

Ro 37 c 2 dt 3j(Ro) Jr dt 



37(^0 

V r 
37(^0 



d 

at 



Ro 



£ / pr 2 dr - Rxp{R{)R\ + R oP (R )R 2 



R2 ° a da(Ro) 



° 2R R + -^d t a 



7(^0 



7 2 ( J R, 



-R 1 p(R 1 )R 2 1 +R p(R )Rl 



(24) 



In the above derivation we pulled u r out of the integration assuming that it equals V r , as 
appropriate in the thin shell approximation. By substituting equation ( |l~8| ) into equation ( p4| ) and 
making use of the following two relations 



V r Ro 



V r R 1 



3c 2 
4 7 2 (iio 



(25) 
(26) 



we get 



f Rl pc 2 f dp u r dp 
Jr 4:jp \dr c 2 dt 



4^0 2 
Pic 



3 7 2 (^o 
—a 



3 7 3 (i?o 



- P {Rv)c 2 



Rl 



3 7 2 0Ro) 



V r V T ■ V T - a 



Rl 



3^(R ) 



V r da(R ). (27) 



Thus, equation (|2^) is now changed to 

■R^ 2 



d t V r 



a 



+ 



Rq 



pic[u r (Ri 
Ri 



V r 



-a 



3 7 2 ( J R C 



'Ro 

-pic 2 +a~ l 



{aR 2 )~ l V r [ Hl [V ■ ( 7P u T )]r 2 dr - (0R 2 )- 1 u r [V ■ ( 7 pu T )]r 2 dr 

JRn JRo 



3 7 3 (i? c 



-p(R )c 2 + 



1 



3 7 2 (R 



■v r v T ■ v T 



+- 



1 



-v r da(Ro). 



(28) 



37 3 (^o) 

In order to close the final equations, we can only have one free variable for the radial velocities. 
We use the approximation that u r is constant across region 2 with V r = u r (Ri), as appropriate 
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under the thin shell approximation. Hence, the first two terms in equation (|28| ) both vanish, and 
we end up with the following equation, 

37^(^0) 3t{Ro) 3t{Ro) 3r\Ro) 

Since 

dq(Ro) V r 

equation (|9|) can be rewritten as 

d a (R ) = -2a- 1 1 (R Q )p 1 c + \a- l p{R Q )c + V r V T ■ V r . (31) 

2 2c 

Because p/p 4 ' 3 is a constant across region 2, we obtain the following relation 

^o)=^.2^^)g. (32) 



Using this result, we can rewrite equation (|3l|) as 

3 3 / 4 _ 

2 l/2 a 7 l/2( i?0 ) c l/2 1 2C 



c\ 7 (# ) = ~2a 1 (R )p lC +—-a + — — V r V T • V r . (33) 



For the tangential velocity we have 

$V T (M) = -^V r -2^V r + CT - 1 ^ 2 7 (i? 1 )p( J R 1 )u T ( J R 1 )[ J R 1 -n r ( J R 1 )] 



+ C j- 1 7(i?o)p(i?o)uT(^o)K(fio) - i?o] - (vRlr 1 ~/p^^r 2 dr 

JRo r 



R 47p V C 2 dt , 

Ri rRi 



(aRl)- 1 f 1 7p(u T • Vu T )r 2 dr - (ai^)" 1 / * u T V • (jpu T )r 2 dr. (34) 



Apparently the last two terms in the above equation are nonlinear and can be ignored. By 
substituting equations (|l3|) and ( fl8| ) into the above equation, we get 

Rl u r u T 2 
7/0 r ar 



d,y T = ir-^^-j Pi c\u t (r 1 )-v t } + v t v t -Vt-(<'RI)- 1 J"' 



-(^r'T^fv^+^t)^, (35) 



R 4 7p V c 2 3t 

The second term on the right-hand-side of equation ( |35| ) is nonlinear and can be neglected. In 
the thin shell approximation, the third term on the right-hand side of equation (|35|) can be 
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approximated as — VtVt / Ro- Using the shock jump conditions at shock 1 and making use of the 
fact that the tangential velocities must be continuous across the shock front, we obtain 



u T {R x ) = -UriR^iyrRi). 

Based on these considerations, equation ( |35| ) can be rewritten as 

2 



(36) 



dtV T 



'- 1 (f )' 



Vr 



Pi c[-MRi)(VtRi) ~ V r ] - -f V T 



2\-l 



R i pc 2 



r dr. 



Next we consider the integration term in the above equation which includes 
f Rl pc 2 f Rl c 2 



(37) 



rill r 2 



/ 1 i-2^T{l P )]r 2 dr- f l P(y Tl ) r *dr 

Jro 3i z Jro 3y 

rRi 

V T / ^pr 2 dr-^(R 1 )p(R 1 )R 2 1 (V T Ri 
Jr 



3 7 2 (i?o) _ 

+"f(R )p(R )R 2 (V T Ro) 



c 2 V T7 (i? ) f Rl 



^ypr 2 dr 



3 7 3 (i? ) Jro 
R 2 V T a + 2aR (V T R ) - 1 {R 1 )p{R l )Rl{V T Ri) 

„2 



3 7 2 (i?o) 
+ 1 (R )p(R )R 2 (y T Ro) 



3 7 3 (« 



o 



(38) 



In the last pair of square brackets of equation (|38|), the second term is much smaller than the 
fourth term and so it can be neglected. Another integration term is 



Rl pc 2 ut dp 2 

r 4 7 p c 2 dt 



/ — -~(9f lnp)^pu T r 2 dr 

Jrq 

l r Ri 
— ^— -r[d t kip(Ri)] / -fpu T r 2 dr 
4t{Ro) Jro 

( dtPi 

4 7 2 (^o) V Pi 



aR%V T - 



(39) 



Because dtlnp does not change much across region 2, we took it out of the integration in the 
above derivation. 



By substituting equation ( |3l| ) into equation (p9|), we get 



[ Rl pl_^L d ± r i dr 
I Bo 4 7 p c 2 dt 



( dtp! 



R 2 R 2 i 

- a ° , PicV T + . ° . p{R )cV T + ~. — oTT, , i 
7 2 (i?o) 4 7 3 (i? ) 4 7 2 ( J R ) V Pi 



(tR$V t . (40) 
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Now, by substituting equations (pq) and (E9) into equation 



BtVi 



-a 



Pi 

,,2 



-a 



p(Rq) 

l(Ro) 
-V T cr + 



VtRq — o 1 PicVt 



-Vtj{Rc 



we get 



i?0 
1 



3 7 3 (^o 



(dtPi\ 
4 7 2(i? ) { Pl J 



(41) 



In deriving the above equation, we made the assumption that surface irregularities due to 
variations in the thickness of the shock layer is of higher order than irregularities due to the bulk 
displacement of the shock layer (regions 2 & 3 in Figure 1), so VtRo = Vt-Ri = Vyi?2- This is 
appropriate under the thin shell approximation. Also the last term in the above equation is much 
smaller than the third term for the wind case and is equal to zero for the ISM case, so can be 
neglected. If we now substitute equation (32) into the above equation, we get equation (^) in §2. 

The derivation of the perturbation equations of region 3 is very similar to that of region 2. 
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Fig. 1. — Structure of the forward/reverse shock system. 




Fig. 2. — Real and imaginary parts of the six eigenvalues in the wind case as functions of I with 
7c = 500. 
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Fig. 3. — Real and imaginary parts of the six components of different eigenvectors (modes) in the 
wind case with 7 C = 500. Each row corresponds to one component of the eigenvector. Different 
line types correspond to six different eigenvectors: the thick solid line refers to the eigenvector of 
Ai, the thick dashed line to A2, the thin solid line to A3, the thin dashed line to A4, the dotted line 
to A5 and the dashed-dot line refers to A6- 
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Fig. 4. — Evolution of the perturbation variables in the wind case with 7 C = 500. Four different 
line types correspond to / = 5, / = 50, / = 500 and / = 5 x 10 3 . 
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Fig. 5. — Evolution of the perturbation variables in the ISM case with 7 C = 500. 



